clear all;
clc;
tstart=1; %start time
tmax=1; %max time of simulation
time=tstart:tmax; %time vector

n_iterations=1;
%nparticles = 1000;
nparticles = 10;
npos=2*length(time)+3; %size of pos arrays = 2t+3
pos=zeros(npos,npos);
pos_x=1:npos;
pos_y=1:npos;

nmid_rel=tmax+2; %relative mid position in the below a arrays which corresponds to zero position in the absolute sense
nstart_rel=0;%relative start position in the below a arrays which corresponds to least -ve position in the absolute sense
nend_rel=0;%relative end position in the below a arrays which corresponds to max +ve position in the absolute sense

a_n_u_prev=zeros(npos,npos); %a's for spin up at each N at previous time
a_n_d_prev=zeros(npos,npos);%a's for spin down at each N at previous time
a_n_u_cur=zeros(npos,npos);%a's for spin up at each N at current time
a_n_d_cur=zeros(npos,npos);%a's for spin down at each N at current time

%%%%%%%%%%%%%%%%%%%%%% SIGMA CODE %%%%%%%%%%%%%%%%%%
%sigma=zeros(1,npos); %the standard deviation of the spatial distribution
sigma=zeros(1,length(time));
sigma_up=zeros(1,length(time));
sigma_down=zeros(1,length(time));

avg_pos=zeros(1,length(time));
avg_pos_up=zeros(1,length(time));
avg_pos_down=zeros(1,length(time));

sig=zeros(1,length(time));
sig_up=zeros(1,length(time));
sig_down=zeros(1,length(time));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% END SIGMA %%%%%%%%%%%%%%%%%%


% T=tmax/2:100:tmax;
T=time;
pos_n_counter=0;
pos_m_counter=0;
a=0;
epsilon = 0.025;
p_n_tmax = zeros(npos,npos);
p_n_u_tmax = zeros(npos,npos);
p_n_d_tmax = zeros(npos,npos);

p_n = zeros(npos,npos);
p_n_u= zeros(npos,npos);
p_n_d= zeros(npos,npos);

%%%%%%%%%%%%%% INIT POSITIONS %%%%%%%%%%%%%
for posx_iter=1:npos
    for posy_iter=1:npos
        pos(posx_iter,posy_iter) = sqrt(posx_iter^2 + posy_iter^2);
    end;
end;

%%%%%%%%%%%%%%%%% END INIT POSITIONS %%%%%%%%%%%%%


for iter=1:n_iterations
    display(iter);
    a = 1/(sqrt(2));
%     a = rand;
    a_n_u_prev(nmid_rel,nmid_rel)=a;
%     a_n_d_prev(nmid_rel,nmid_rel)=(sqrt(1-a^2))*i;
    a_n_d_prev(nmid_rel,nmid_rel)=a*i;
    
    
    %for all times
    for t=time
        display(t);
        nstart_rel = nmid_rel - t;
        nend_rel = nmid_rel + t;
        if(t ~= 1)
             a=1/sqrt(2);
%              a=(1/(sqrt(2))) + epsilon*rand;
%             a = epsilon*rand;
        end;

        %for all possible positions for this particular time, fill up the a
        %arrays for the current time
      pos_n_counter=nmid_rel;
%     for pos_n_counter=nstart_rel:nend_rel
            
        for pos_m_counter=nstart_rel:nend_rel
%             display(pos_n_counter);
%                   display(pos_m_counter);
%                 if(rand <= 0.5)                    
%                     a_n_u_cur(pos_n_counter,pos_m_counter) =  ( sqrt(1-a*a) * a_n_d_prev(pos_n_counter-1,pos_m_counter) +  a * a_n_u_prev(pos_n_counter-1,pos_m_counter) );
%                     a_n_d_cur(pos_n_counter,pos_m_counter) =  ( sqrt(1-a*a) * a_n_u_prev(pos_n_counter+1,pos_m_counter) -  a * a_n_d_prev(pos_n_counter+1,pos_m_counter) );                    
                    
%                 else
                    a_n_u_cur(pos_n_counter,pos_m_counter) =  ( sqrt(1-a*a) * a_n_d_prev(pos_n_counter,pos_m_counter-1) +  a * a_n_u_prev(pos_n_counter,pos_m_counter-1) );
                    a_n_d_cur(pos_n_counter,pos_m_counter) =  ( sqrt(1-a*a) * a_n_u_prev(pos_n_counter,pos_m_counter+1) -  a * a_n_d_prev(pos_n_counter,pos_m_counter+1) );                    
% 
%                 end;

                %record the probability densities of up,down and all particles
%                 p_n_u_tmax(pos_n_counter,pos_m_counter) = (abs(a_n_u_cur(pos_n_counter,pos_m_counter)))^2;
%                 p_n_d_tmax(pos_n_counter,pos_m_counter) = (abs(a_n_d_cur(pos_n_counter,pos_m_counter)))^2; 
%                 p_n_tmax(pos_n_counter,pos_m_counter) = (abs(a_n_u_cur(pos_n_counter,pos_m_counter)))^2 + (abs(a_n_d_cur(pos_n_counter,pos_m_counter)))^2; 
            end;
%         end;
        
        %record the probability densities of up,down and all particles
        p_n_u_tmax = abs(a_n_u_cur).^2;
        p_n_d_tmax = abs(a_n_d_cur).^2; 
        p_n_tmax =  p_n_u_tmax + p_n_d_tmax;     
        
        p_n = p_n + p_n_tmax;
        p_n_u = p_n_u + p_n_u_tmax;
        p_n_d = p_n_d + p_n_d_tmax;

        
        %%%%%%%%%%%%%%%%%%%%%%%%%%% SIGMA CODE %%%%%%%%%%%%%%%%%%%%%%
%         if(iter == n_iterations)
%             sigma(t) = (sum( sum((pos-nmid_rel).^2 .* p_n_tmax)) - sum(sum(pos-nmid_rel.*p_n_tmax)).^2 );
%             sigma_up(t) = (sum( sum((pos-nmid_rel).^2 .* p_n_u_tmax)) - sum(sum(pos-nmid_rel.*p_n_u_tmax)).^2 );
%             sigma_down(t) = (sum( sum((pos-nmid_rel).^2 .* p_n_d_tmax)) - sum(sum(pos-nmid_rel.*p_n_d_tmax)).^2);
% 
% 
% %                 avg_pos(t) = sum((pos-nmid_rel).* p_n_tmax);
% %                 avg_pos_up(t) = sum((pos-nmid_rel).* p_n_u_tmax);
% %                 avg_pos_down(t) = sum((pos-nmid_rel).* p_n_d_tmax);
% 
%             sig = sig + sigma;
%             sig_up = sig_up + sigma_up;
%             sig_down = sig_down + sigma_down;
%         end;
        %%%%%%%%%%%%%%%%%%%%% END SIGMA CODE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
     
        
        %set the a arrays of prev time with values of that of the current
        %time and a arrays of current time with zero values
        if(t <tmax)
            a_n_u_prev = a_n_u_cur;
            a_n_d_prev = a_n_d_cur;
            a_n_u_cur = zeros(npos,npos);
            a_n_d_cur = zeros(npos,npos);
        end;
        
        if(iter ~= n_iterations) 
            p_n_tmax = zeros(npos,npos);
            p_n_u_tmax = zeros(npos,npos);
            p_n_d_tmax = zeros(npos,npos);
        end;
        
    end;
   
end;

p_n = p_n/(n_iterations*tmax);
p_n_u = p_n_u/(n_iterations*tmax);
p_n_d = p_n_d/(n_iterations*tmax);
%

% sig = sig/(n_iterations*tmax);
% sig_up = sig_up/(n_iterations*tmax);
% sig_down = sig_down/(n_iterations*tmax);


% figure(9)
% plot(T,sig,T,sig_up,T,sig_down);
% xlabel('T:1000');
% ylabel('\sigma^2(t) for charge and up-spin,down-spin');
% title('\sigma^2(t) for charge and up-spin,down-spin, for \epsilon = 0.3 and a=1/sqrt(5)');
% grid on;


figure(10)
surf(pos_x-nmid_rel,pos_y-nmid_rel,p_n);
title('plot of P_n(t) for 2D Quantum hadamard random walk');
xlabel('X axis');
ylabel('Y axis');
zlabel('P_n(t)');
grid on;

% figure(11)
% surf(pos_x-nmid_rel,pos_y-nmid_rel,p_n_u);
% title('plot of P_n(t) for 2D Quantum hadamard random walk - up-spin');
% grid on;
% % 
% figure(12)
% surf(pos_x-nmid_rel,pos_y-nmid_rel,p_n_d);
% title('plot of P_n(t) for 2D Quantum hadamard random walk - down-spin');
% grid on;